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Abstract 

The Schrodinger-Langevin equation (SLE) is considered as an effective open quantum 
system formalism suitable for phenomenological applications involving a quantum sub¬ 
system interacting with a thermal bath. We focus on two open issues relative to its 
solutions: the stationarity of the excited states of the non-interacting subsystem when 
one considers the dissipation only and the thermal relaxation toward asymptotic distri¬ 
butions with the additional stochastic term. We first show that a proper application of 
the Madelung/polar transformation of the wave function leads to a non zero damping 
of the excited states of the quantum subsystem. We then study analytically and nu¬ 
merically the SLE ability to bring a quantum subsystem to the thermal equilibrium of 
statistical mechanics. To do so, concepts about statistical mixed states and quantum 
noises are discussed and a detailed analysis is carried with two kinds of noise and po¬ 
tential. We show that within our assumptions the use of the SLE as an effective open 
quantum system formalism is possible and discuss some of its limitations. 

Keywords: open quantum system, Schrodinger-Langevin equation, thermal 
relaxation, stationarity 


1. Introduction 

In classical mechanics, the influence of a thermal environment (bath) on a Brownian 
particle (subsystem) is well described by the Langevin dynamics within the Newtonian 
framework. The subsystem thermalization is obtained from the balance of two forces 
(friction and stochastic) which generate irreversible energy exchanges between the two 
systems. To search for the corresponding description in quantum mechanics is a crucial 
issue both for the understanding of quantum fundamentals and in many branches of 
applied physics (where the quantum systems can never be isolated), such as in quantum 
diffusion and transport mmm, quantum optics laE], heavy ion scattering PEllH], 
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quantum computers and devices pnoKii]. Unfortunately, the Langevin dynamics - or 
more generally energy dissipation - cannot be introduced easily in the common quantum 
formalism, as no direct canonical quantization of a Hamiltonian can describe irreversible 
phenomena m- 


In the present work, we focus on a possible Langevin-like extension of the fundamen¬ 
tal Schrodinger equation, the so-called Schrodinger-Langevin equation (“SLE”) 


ih 


d'ijj{x, t) 
dt 


Hq + hA(^S{x,t) — [ ij}*S{ xA)' 4’ dx^ — xFR{t) 




( 1 ) 


where A is the friction coefficient, S the (real) phase of the wave function, a 

fluctuation operator and F[q the usual isolated-subsystem Hamiltonian, 

Ho = + Uext(x). (2) 

The SLE was first proposed by Kostin m from an identification with the Langevin 
equation for Heisenberg operators, the so-called Heisenberg-Langevin equation (“HLE”), 

P = Fe^tiX)-AP + FR{t) and X = Plm. (3) 

The latter is derived within the common subsystem plus bath approactQ from a simple 
model of the bath [EiiigiiT! - a thermal ensemble of oscillators linearly coupled to the 
subsystem - and has proven to be a suitable framework to study Brownian motion. The 
practical application of the HLE is nevertheless limited by its non-commutating operator 
nature. The SLE has also been derived within many other frameworks and non-standard 
quantization procedures to describe either pure dissipation [ig na [23 EH i22i [231 [21] 
or the case of a Brownian motion [SHESIETI. The SLE includes a thermal fluctuation 
term —xFji(t) and a dissipative term under its hydrodynamic formulation |20[ [2H [25] 

hA{S{x,t)-{S)), (4) 


where the phase S is chosen according to a prescription which will be discussed in 
Sec. Even though the dissipative term is nonlinearly (logarithmically) dependent on 
the wave function, it still corresponds to a linear ohmic friction (i.e. proportional to the 
particle velocity). A nonlinear friction can be obtained by extending Kostin derivation 
to a nonlinear coupling |28[ 129) . The SLE exhibits interesting properties: unitarity is 
preserved at all times [21], the uncertainty principle is always satisfiec^ [HI EH [33] 
and the superposition principle is violated due to the nonlinearities (which might not 
be a problem per se for dissipative equations [31E5]). Thanks to its straightforward 


^In the common approach, the subsystem plus bath is considered as a whole conservative system. By 
integrating out the bath degrees of freedom, one obtains the dissipative evolution of the subsystem m- 
^As opposed to other models like the Caldirola-Kanai equation [^ without fluctuations, the Wigner- 
Moyal equation with classical Fokker-Planck terms or the quasiclassical HLE |14| . 
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formulation - in principle only two “classical” parameters need to be known: the friction 
coefficient A and the bath temperature Tbath ^ and its numerical simplicity, the SLE can 
be considered as a solid candidate for effective description of open quantum systems 
hardly accessible to quantum master equations [HEH]. Indeed, in a number of complex 
applications, defining the bath/interaction Hamiltonian and calculating the Lindblad 
operators without too many approximations is rather complicated, and some effective 
approaches - possibly of the Langevin type - are unavoidable naEZj. Because there is no 
established connection between the SLE and the standard quantum master equations, the 
SLE is different from the stochastic Schrodinger equation (SSE) developed to mock the 
evolutions given by the quantum master equations [38l [39l SO] ■ As stochastic equations 
based on pure state evolutions, they nevertheless share the same philosophy to perform 
an average over a large ensemble of initially identical subsystems to recover the statistical 
mixed state describing the subsysterr0 The expectation value of an observable operator 
O is then given by 

-j ^stat 

((tp(t)ldlip(t)}\ = lim - , (5) 

\ / stat ^^stat^oo TT-stat 

r=l 

where the pure state is given by the r*^ realization of the stochastic evolu¬ 

tion. The numerical costs remain quite reasonable in comparison to the common density 
matrix approach which is highly expensive when the Hilbert space associated to the 
subsystem is large. 

Before considering any actual application to phenomenology, some questions and 
issues remain to be explored about the solutions of the SLE and its thermal relaxation. 

The study of its solutions without stochastic term has been carried out in many 
specific cases, either analytically na ED ESI mi Ha aa im sa He] or numerically [241 E3l 
mmaiisiiiaEQiED. Along these analysis, it has been advocated that the stationary 
eigenstates of Hq are also stationary states of the equation |18l I43j . This behavior is in 
contradiction with what is expected from damped quantum systems da Ea Eg ED El!- 
As an answer to this expectation, we will first show in Sec.j^how a proper application of 
the Madelung/polar transformation of the wave function results in the damping of these 
states. The purely dissipative SLE has already been applied in quantum chemistry [23] 
and heavy ion scattering jSa E3] • 

However very few studies including the stochastic term have been carried out and 
these are moreover mere comparisons with the solutions of the HLE. Kostin m first 
observed that for a free particle plane wave, the SLE and HLE lead to the same solu¬ 
tion. Then, Messer [53] studied the evolution of a Gaussian wave packet in the free and 
harmonic cases. In the free case, he showed that the evolution differs from the HLE so¬ 
lution, highlighting that the SLE and HLE are not strictly equivalent. In his calculation. 


^Because of the statistical nature of the bath-subsystem interactions, the subsystem must be described 
by a mixed state, which includes not only probabilistic information about the observable measurements 
but also about the state itself. The common tool to describe a mixed state is the density matrix operator. 
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Messer used a white noise for the stochastic force - which is questionable - and 

assumed that the SLE naturally leads to the thermal equilibrium predicted by statistical 
mechanic^ (the Gibbs state). The latter is characterized by a Boltzmann distribution 
of the uncoupled subsystem energy states {En}n=o,i...- 

( ^ 

Pn OC exp —- , 

V ^^bath / 

where k is the Boltzmann constant and pn the population (or “weight”) of the eigen¬ 
state of energy En, and is generally expected at the weak coupling limit for a quantum 
subsystem in interaction with a heat bath [SZIIMIISZI. The so-called weak coupling limit 
or Brownian hierarchy is achieved when the relaxation time of the subsystem (~ 1/^) is 
much larger than its natural period of oscillation and than the typical correlation time 
a of the microscopic interactions between the bath components and the subsystem. To 
our knowledge it has never been proven or tested that the SLE actually admits such an 
asymptotic distribution in a dynamical manner. The main contribution of the present 
work is precisely to study the thermal relaxation given by the SLE with different ID 
subsystems using either a white or a colored noise and thus to test Messer’s assumption 
as a byproduct. To this end, in section we discuss different possible choices for the 
stochastic force Fji{t) (which will be assumed to be a c-number). In sectionwe study 
the equilibration given by the SLE with a ID harmonic potential I4xt = K 2r/2. To do 
so, we first show analytically that the SLE brings a Gaussian wave packet to the Gibbs 
state if one uses a specific white noise. We generalize this result to other initial states 
through the numerical resolution of the SLE with the Grank-Nicolson scheme. Similarly, 
we then explore the equilibration resulting from a colored noise and demonstrate that 
it leads to the correct equilibration in the weak coupling limit. In section we extend 
the numerical simulations to a linear potential I4xt = Fi \x\/2 to study the equilibration 
given by the SLE within a non-harmonic situation. 

Though the thermal equilibrium of statistical mechanics is only expected at the weak 
coupling hmi10 (^4 <C {cuo, l/^'})) where wq is the characteristic frequency associated to 
the potential I4xt0 the intermediate {A < {wqjI/o'}) and strongly coupled regimes 
[A > {cuo, l/c}) are also investigated. 

Although we perform the numerical simulations with the dimensionless SLE for sim¬ 
plicity - i.e. with natural units h = m = K = Ki=ujQ = k = l and dimensionless 
variables x, t,. . .[^ - the analysis shows that the important dimensionless ratios gov- 


( 6 ) 


"^Within these assumptions, the SLE has already been applied to atomic diffusion in solids m- 

® Generally, a quantum subsystem in interaction with a heat bath is expected to reach a thermal 
equilibrium, where the components of its energy spectrum are shifted and broadened |,S7[ 1561157 | . These 
spectrum modifications become negligible at the weak coupling limit and one expects the thermal equi¬ 
librium of statistical mechanics. Note also that the equations of evolution, such as the quantum master 
equations and the HLE, are derived under the assumption of a weak coupling. 

^ujo = for the harmonic potential and u>o = ^ for the linear potential. 

^The dimensioned values of x, t, A, Fr and Hq can be obtained by multiplying our dimensionless 
values respectively by y/himujo, I/oiq, ujq, \/mhujg and hujQ, where uq is specific to each type of potential. 
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erning the physics are Ajuj^ and feTbath/(^<^o) for both potentials. In the main part 
of the text, it is then assumed, for the purpose of compactness, that A and Tbath are 
respectively measured in units of wq and huj^/k - and thus correspond to these ratios 
while times are measured in units of and energies in units of Hojq. For all analyti¬ 
cal calculations, we stick however to the International System of Units, for the sake of 
clarity. 

Along this work, we will extract some effective temperature Tguh reached by the 
subsystem from the asymptotic weights {pn}n=o,i..., mostly by fitting a Boltzmann dis¬ 
tribution oc exp (—En/kTsuh) to this equilibrium distribution. In some specific situations 
we will observe that Tsuh differs from the bath temperature Tbath (which is concretely 
defined as the temperature entering the noise correlation); it is precisely one of the main 
goals of our study to identify under which conditions such a departure happens. 

2. A well defined prescription for the friction term to obtain eigenstates 

damping 

Though theoretically the fluctuation and dissipation aspects cannot be dissociated, 
it appears that in some specific studies only the damping is considered [HI [511 [52l 
123] • Unfortunately, the dissipative part of the SLE 0 suffers from ambiguities and 
some prescription is required to obtain an unambiguous definition. Its main non-linear 
ingredient is the real phase S{x,t), defined by the wave function decomposition 

(7) 

where R{x,t) is the real amplitude. S{x,t) is indeterminate at the wave function nodes 
and, in general, multivalued (defined modulo 27r). 

In the literature jlHl (HI SSI EE], a common prescription adopted for real il) is to 
require the phase S{'ip) to be zero (and thus continuous at the nodes of ip) while R{ip) 
is taken as a real - positive or negative - function. This prescription has led to the 
conclusion that the stationary eigenstates of Hq are also stationary states of the SLE, 
as the dissipation term identically vanishes. For the sake of describing time-dependent 
situations, a corresponding prescription has however to be adopted for any complex 
Ip as well. It is easily seen that such an analytical continuation unavoidably has one 
branch cut in each half complex-plane, both of them starting from the origin. Taking 
for instance those branch cuts along the imaginary axis corresponds to 

S{ip) = arcta,ii{'is{ip)/ilt{ip)). (8) 

with finite values of the friction potential in the SLE and then finite damping. Therefore, 
an infinitesimal modification of ip (associated to a slight deviation from real axis to 
complex plane) leads to a large variation of the associated damping of the quantum 
state, which is the sign of an ill-defined model. 

We propose to use instead the “polar” or “Madelung” prescription, where one defines 
R{x, t) in eq. Q as the module of the wave function, i.e. a real positive function. In 
practice, one could use the local argument of the wave function ATg(ip) - a well-defined 
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function with a branch cut along the negative real axi^- to determine S{x, t). However, 
the limitation of the Arg values to a ] — vr, +7r] interval (as illustrated in Fig. would 
lead to discontinuities of the dissipative term with unphysical effect^ To avoid these, 
we suggest to build the phase S{x, t) on a spacial grid of step dx following the recursive 
law 

S{x + dx) = S{x) + dS{x) where dS{x) = Aig['ij;{x + dx)/'ip{x)], (9) 

starting from an arbitrary space point of reference “0”. This leads to 

j 

S{j X dx) = S{0) + dS{k X dx). (10) 

k=l 


The chosen value of the multivalued 5'(0) is of no importance thanks to the regulator 
— (S) and can therefore be taken to Arg(?/;(0)). 

S(x)-<S) lAaW and its S(x)-(S) 




Figure 1: The dissipative term S{x) — {S) 
corresponding to the plane wave ipi^) = 
obtained with the argument (Arg) function 
(dashed line) and with the recursive method 
(solid line). 


Figure 2: The real '03 harmonic eigenstate 
(solid line, magnified for the plot sake) and 
its corresponding phase S{x) — {S) (dashed 
line) with the recursive method in the polar 
TT prescription. 


The polar prescription implies singular phase shifts +Tr at the wave function nodes 
as shown for instance in Fig. Not only are these discontinuities theoretically allowed 
(thanks to the phase indeterminacy at the nodes), but they also have a convenient 
physical consequence: the eigenstates of Hq are not stationary states of the SLE anymore. 
Indeed, for the excited eigenstates {ipnjn^i of Hq the friction term becomes a step 
potential which generates correlations between eigenstates and results in damping. To 
show the latter assertions, let us assume that an initial wave function '0(0) is equal to 
an eigenstate 0m>i, i-e. 0(0) = ^ Cn(O)0n with Cn(0) = dnm- The SLE without the 


®In many programming languages, such a choice of the argument of a complex number 2 is given by 
the function atan2, i.e. Arg( 2 ) = atan2(Q'(2), 5R(2)). 

®The invariance under the multiplication of the wave function by a simple phase factor would be 
broken. 
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stochastic term yields, 


Cn = -■^(V'nli^olV') - (V’nK-S* - (S'))|V') 

h 

= -^EnCn-iA'^Ck [{S - {S))'ll)*^'4>kdx . (11) 

For symmetric potentials for instance, one can show that if 'il’k=m has an odd (even) 
parity, then the integral is finite and thus the transition Cm —?• Cn is allowed at very 
small times for all ipn with even (odd) parities. Moreover, the smaller the difference 
|n — m|, the larger the transition rate, which is consistent with the typical behavior of 
transition matrix elements entering e.g. the Fermi Golden Rule. Last but not least, the 
transition rate ton = m—1 is larger than to n = m+1, which is consistent with damping. 
At larger times, these transitions and the damping can be observed numerically (see for 
instance Fig. [^. 



Figure 3: Evolution of the eigenstate weights Pn=o,i ,2 = l(V'n|V'(^))P the dimensionless SLE 
without fluctuating term from an initial first excited state, with the “polar” prescription (solid 
lines) and the “arctan” prescription (dashed lines). 

Both the “arctan” and “polar” prescriptions are mathematically correct and the 
choice between them should be physically motivated. The stationarity of the Hq eigen¬ 
states in the corresponding dissipative situation remains an open question within the 
open quantum system framework [59] . As illustrated in Fig. the SLE can reproduce 
both situations, thanks to the two prescriptions presented in this section. From the 
perspective of the common quantum master equation ISZI, the Lamb shifted energy lev¬ 
els acquire hnite lifetimes (finite widths), implying that the polar prescription is better 
suited for robust phenomenological studies. Let us hnally stress that the choice of the 
prescription is of little importance when the fluctuations are considered, as they drive 
the state away from any given eigenstate. 
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3. The possible noises for the stochastic term 

We now focus on the stochastic aspect of the SLE Q. Similarly to other Langevin- 
like equations, the noise term simulates the many collisions (or couplings) that 

the subsystem undergoes with the particles of the bath. It is generally taken as a 
homogeneous Gaussian random process, independent of the subsystem position, and 
described by its average and covariance function. The random direction of the many 
collisions always yields a zero average. The covariance function C is usually given by 
a fluctuation-dissipation relation. The latter relates the noise covariance, the friction 
coefficient A and the bath temperature Tbath in order to obtain a balance between the 
fluctuation and dissipation aspects. The balance is correct if the subsystem distribution 
at equilibrium is Boltzmannian oc exp (—E„/Tsub) and if the temperature reached by the 
subsystem T^nh is equal to Tbath- As the fluctuation-dissipation relation corresponding 
to the SLE has never been determined to our knowledge, our method will consist in 
borrowing such relation from another framework (discussion in this section) and then 
evaluating its consequences for the SLE (work done in next sections). As the SLE might 
be the counterpart of the HLE in the Schrodinger representation, the covariance func¬ 
tions derived within the HLE framework might be suited. Within this framework, the 
noise operator is built from the initial bath position and momentum operators whose 
non-commutative property leads to the main differences with the classical case. 


Senitzky |15) first proposed an HLE - for a general bath linearly acting on a harmonic 
subsystem (with natural frequency cjq) “ where the noise operator is described by a white 
noise autocorrelation, 


hujQ _(iWO_le/ N 

2 exp(hwo/A:Tbath) - ij 

( 12 ) 

where 6 is the Dirac distribution. As in the classical case, a white noise covariance 
means that there is no correlation between the collisions, which results in a Markovian 
process. This covariance has been used by Messer [5lj in its analytic comparison of the 
HLE and SLE solutions. The first term of the RHS bracket in the definition of Csen 
corresponds to the zero point fluctuations of the subsystem. This term is required within 
the HLE framework for the canonical commutations to hold at Tbath = 0, as shown by 
equation (52) in |15] . However, within the SL framework, the zero point fluctuations 
appear naturally in the wave function, so that they do not need to be included in the 
noise operator for the canonical commutations to hold. Therefore, this term becomes 
unnecessary and the white quantum noise autocorrelation Gwhite will be defined as 


(^FR{t)FR{t + r)^ = Csen(r) with Csen{r) := 2mA 


C'white('’‘) := B 5{t) with B := 2mA Eq 


coth 


To 


/cT 


bath 


-1 


(13) 


where Eq = hujQj2, the zero point energy. In Sec. we will show that the fluctuation- 
dissipation relation (13) indeed allows to reach asymptotically a thermal distribution of 


states when one uses a white noise and a harmonic potential. 
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However, Li et al. [60] pointed out an important weakness in the derivation of (12) 
They also claimed that the colored quantum noise 


FR{t)FR{t + t)) = 


m 


TT 


f 


huj coth ( ——-] cos(ujt) + i sin(a;r) A duj. 

\ 2A^7bath / 


(14) 


first derived by Ford et al. [T6| , is the only one able to drive a general subsystem to the 
correct thermal equilibrium via the HLE. For now, the latter assertion has only been 
demonstrated in a limited form for harmonic and nearly harmonic oscillators [561IHT] . 

Actually, in order to get rid of the contribution from the bath zero point fluctuations 
- which was first judged physically unjustihed -, Ford et al. first derived a colored 
quantum noise under the form of the normal product 


N[FR{t)FR{t + r)]) = C'coiored(T) := 


hui 


2mA 


TT 7o exp(;iw/A;Tbath) - 1 


cos(wr) doj. (15) 


Both colored noises lead to a non-Markovian process even if the friction is memory-less. 
As pointed out by Gardiner |56|, the correct choice of spectrum depends on what is 
actually measured to find it: e.g. in absorption measurements one gets the black body 
radiation Planck spectrum corresponding to (15), whereas in Josephson junction noise 
current measurements 


one gets the linearly rising spectrum at high frequencies cor¬ 
responding to (HI- Within the SSE framework, the choice between the different noises 
is also intensively discussed, for instance to obtain the correct thermal equilibrium of 
a non-Markovian master equation m or the correct positivity property for the Bloch- 
Redfield master equation [63] . 


In general, the practical application of the HLE is limited by its non-commutating 
operator nature. Although questionable |64j . a common approximation [651 IMlETllMl 
ESlIZQlCI] is to abandon its operator character and to replace the non-commutating 
q-number noise by a c-number noise while taking the same power spectrum. One then 
obtains a quasiclassical Langevin equation which leads to a reasonable description for 
systems which are nearly harmonic and to possible violations of the Heisenberg rela¬ 
tions. Within the SLE, we only need to assume that the noise operator can be taken 
as a commutating c-number. The latter assumption was actually already implied in 
Kostin’s derivation of the SL random potential m and does not lead to a violation of 
the Heisenberg relations [3I1E21E3]. 

In the present paper, we adopt the same assumption: Fr will be considered as 
Gaussian stochastic c-numbers of zero average while the autocorrelation C will be taken 
according either to dehnition (13) (white noise) or definition (colored noise). We will 
then test the ability of the SLE to bring a subsystem to thermal equilibrium. Whereas 


the white noise (13) leads to uncorrelated stochastic forces, the colored noise (15) con¬ 


tains a strong temperature dependence of their correlation time. The latter becomes 
large at low temperatures (oc 1/Tbath) and the Brownian hierarchy/weak coupling limit 
could be broken in such a situation, when A > Tbath (see introduction). In Appendix 


[A] we describe the algorithm used to build these correlated forces numerically. The 
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colored noise (14) would have also been pertinent because of its direct connection with 
the quantum dissipation-fluctuation theorem, but it leads to additional complications 
(imaginary part) and model dependences. Indeed, the noise correlation (14) strongly 


depends on the value of its high frequency cut-off - which evaluation is specific to each 
system - and on the choice of the cut-off shape (Lorentzian, exponential, sharp...). 


4. Equilibration with a harmonic oscillator 

In this section, we study the thermal relaxation given by the SLE (|^ with the 
harmonic ID potential I4xt = mcuQX^/2 and with the white (13) or colored (15) 


noise. 


The harmonic oscillator is a well-known basis to study the properties of an open quantum 
system formalism 


4-1- Wave function evolution during one stochastic realization 

We first focus on the evolution of the wave function during one stochastic realization 
of the SLE. Analytically, one can study the evolution of a general Gaussian wave packet 
Ansatz 


lfA{x,t) = gS ("W(a;-2:cl(i))^+Pcl(i)(a;-a:cl(i))+7(i)) ^ 


(16) 


where a{t) is a complex number related to the wave packet width (Im(Q!) > 0 at all 
times), 7 (t) a complex phase whose imaginary part plays the role of a normalizing 
factor, and Xc\ and pc\ are centroids (and expectation values) in direct and momentum 


space. As shown in Appendix B, inserting (16) in the SLE leads to ordinary differential 


equations for a, aich Pel and 7 , including: 


a + A3f?(a) -|- —a^ -|- = 0 

m 2 


(17) 


and 


Pci 


Pel = -mujQ Xci - Apci + Er , Xci = — 


m 


(18) 


In [Appendix B it is also demonstrated that the solution of eq. ( |17[ ) starting from any 
Q{a(t = 0 )) > 0 tends asymptotically to a{t —)• 00 ) = imuJof2, which corresponds to 
the width of the ground state a = -s/h/mujQ . After some initial relaxation, the general 
solution (16) from any initial Gaussian state is thus the ground state displaced in direct 


and momentum space with a stochastic trajectory obeying the classical equations of 


motion (18). 


This result can be generalized to other initial states through a numerical resolution 
of the dimensionless SLE with the Crank-Nicolson scheme, which requires solving the 
equation 


1 + 


iH{t)At 


ijj{t + At) = 1 - 


iH{t)At 


i^{t) 


(19) 
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l</'(t=0)r 


i</'(t=i)r 


|iA(t=4)p 



|i^(t=10)P 



Figure 4: Typical wave function shape/module evolution toward the ground state shape/module 
during one noise realization. 



Figure 5: Typical wave function phase evolution toward linearity (where the wave function takes 
significant values) during one noise realization. 


with respect to ^p{t + At) at each time step At. For this purpose, the wave function 
has first been discretized on a spatial grid, with the grid size chosen large enough to 
avoid spurious reflections and the grid spacing much smaller then the typical inverse 
wave number. Next, the Thomas algorithm was applied for the inversion of equation 
(19). Note that the stochastic force Tr was considered to be constant on each time step. 
This leads to an effective autocorrelation time ~ At in the white noise case, which has 
however shown to have no practical consequence. Finally, note that the numerical cost 
when considering ensemble averaged observables is proportional to the space-time grid 
and to the number of realizations, i.e. to ngpace x ^time x Ustat (where typically ngpace is 
of the order of the hundreds, ntime and Ustat of the thousands). 


As illustrated in Fig. and we first confirm from observations that after some 
initial relaxation, the general solution from any initial state is the ground state displaced 
in direct and momentum space with a stochastic trajectory. Indeed, a common wave 
function evolution pattern emerges during a noise realization. First, the shape of the 
wave function evolves toward the ground state shape. In parallel, if one starts from an 
initial excited eigenstate, the phase “breaks” at the nodes and evolves toward a linear 
phase in the region where the wave function takes non negligible values (see at t = 17 in 
Fig. for instance). In parallel and until the end of the evolution, the centroid oscillates 
around the potential minimum following a stochastic trajectory along the space axis. 
Some discrepancies to this pattern, coming from numerical instabilities, appear when 
A <C Tbath and when Tbath > 1- 
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4-2. Evolutions of mixed state observables: energy and populations 

We now focus on the evolutions of mixed state observables as given by the statistical 
relation ([^. To do so, we perform the numerical simulations 
^bath = 1 and with the white noise (13) for illustration. 
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in a bath at temperature 


(Ho) 



Figure 6: Solid curves'. Numerical average energy {{Hq)) evolutions for different values of the 
friction coefficient A. Dashed horizontal line. Corresponding HLE theoretical asymptotic value 
given by the exact relation (211. Dashed curves'. Corresponding theoretical evolutions given by 
relation (20). 


We first study the evolution of the average energy which will be written 

{Hq) for simplification, starting from the initial ground state i/jq. Three average energy 
evolutions with friction coefficients corresponding to weak [A = 0.1), intermediate {A = 
0.5) and strong {A = 1) couplings (see introduction) are shown in Fig. 

The average energy evolution rate predicted by Senitzky m within the HLE frame¬ 
work, 

{Ho){t) = + {Ho){t ^ oo)(l - , (20) 

fits well our numerical evolution in the weak coupling case (where Senitzky’s HLE ac¬ 
tually applies) as shown in Fig. Furthermore, the theoretical asymptotic value for a 
quantum harmonic oscillator in thermal equilibrium is given by. 


{Ho){t oo) 


Eq coth 



( 21 ) 


and corresponds to our value {HQ){t — >■ oo) ~ 1.07 when Tbath — 1 for any coupling. 

The second interesting observable is the distribution of the eigenstate weights (pop¬ 
ulations) Pn{t), as given by relation ([^ with the projection operator On = |'0n)('0n|- As 
shown in Fig. their evolutions during the transient phase follow the general expecta¬ 
tion of open quantum systems: the main transitions occur between neighboring energy 


“The grid steps are taken as Aa; = 0.1 and At = 0.01. The covariance parameter of the white noise 
is taken to ct = 0.03; see remark after equation (A.6l. 
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Figure 7: Evolutions of the eigenstate weights Pn=o,...&it) from the initial ground state (left) and 
2"'^ excited state (right) for a friction coefficient corresponding to an intermediate coupling. 


levels. Moreover, they lead to a reshuffling of the weights, such as pn > Pn+i, reached 
after a lapse of time proportional to the relaxation time XjA. 

4--3. Equilibration with a white noise 

Here we show that the equilibrium distribution of the eigenstate weights is a Boltz¬ 
mann distribution oc with = ^bath provided that the fluctuation-dissipa¬ 


tion relation (13) is satisfied if one uses a white noise. In other terms, the subsystem 


equilibrates with the bath if (13) is satisfied when one uses a white noise. We first focus 


on the analytic solution with initial Ganssian wave packets and then show throngh nn- 
merical simulations that these resnlts are universal, i.e. independent of the chosen initial 
state, friction coefficient and temperature. 

4 . 3 . 1 . With Gaussian initial wave packets 


We showed in Sec. 4J. that any Gaussian initial wave packet reduces asymptotically 
to a coherent state with a width a = \ corresponding to the ground state, i.e. 

y TflbJQ 


'ip{t » H oc e ^ 




+i- 


( 22 ) 


We wonld like to evaluate the weight of the different Hq eigenstates il^n in this asymptotic 


^l>. For this purpose, we first reformulate (22) as 


.Ki/gT 


V’ oc e 2 2 

^cl I 

where we have set ^ f and p = —— ^ ^ . Using the identity 


(23) 


+ 00 




n=0 


^ri\ 


(24) 
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where {-f^n}n=o,i... are the Hermite polynomials, as well as the expression for the eigen¬ 
states 




HnjQe-^ 


(25) 


yields 


^ oc e 1 2a j I 2R 1 + 2h \ - 

vn! 


n=0 


(26) 


We thus deduce that the eigenstate weight Pn(a^chPci) for a given realization of the 
stochastic noise is given by 


Pn(Xcl,Pcl) oc 


in I,, 12n 


n! 


_ 1 
e 2 


n! 

(27) 


and one has exactly ^pn = 1- In Sec. 4.1 we have shown that the position Xd and 


momentum pd centroids satisfy the classical stochastic equation of motion (18). When 


the stochastic force autocorrelation is of the form + t)) = B 5{t) (white 

noise) it is known that the distribution of the trajectories (xcb Pci) is 


2 2 2 
, PrI 

—^21 -^■’77^ 


W{xd,Pc\)cce 


(28) 


where Tefj := 2 FA- eigenstate weight, averaged over the fluctuations, will then be 
given by 

Pn= W{xd,Pcl)Pn{XchPcl)dXddpd- (29) 


To determine (29), we use the relation 


Pn(Xcl,Pcl) = ^^^e V J 


After some trivial integration on Xd and pd, one gets that 


ri=l 


(30) 


_ f I (Pcl°)^ ^ ^^0 

W{xd,Pci)e V ^ ^ ydxddpd= 


I hwQ 

d + fcTeff 


(31) 


where the numerators guarantees that for r/ = 0, one has J W{xd,Pci)dxddpd = 1- 
Differentiating n times (31) with respect to p yields 


hwQ / k \ 

^ oc 


Pn = 


1 I hujQ y 
^ fcTeff J 


(32) 
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Following relation (13) and the definition of Tgfj, one finds 


urr ^0 
fcTefr = — 


and 




coth 


1 


huiQ 

2kTy 


bath 


exp - 


- 1 


= exp 


hujo 

^^bath 


-P(2^)—p(-2^) 

The ensuing expression for pn is 

nhuiQ 


Pn oc. e '“^bath (X e *^bath J 


(33) 

(34) 


(35) 


which demonstrates that the distribution of the state weights is Boltzmannian if one uses 
a white noise with autocorrelation (13), with, moreover, T^nh = ^bath for all coupling 
strengths and all bath temperatures. This reasoning can be easily extended to the 3D 
case. 


4-3.2. With other initial states 

We now perform the numerical simulations with the dimensionless SLE to generalize 
the previous results to other initial states. As shown for instance in Fig.[^ the asymptotic 
distribution of the weights is independent of the chosen initial state and perfectly fits 
a Boltzmann distribution. One can determine the temperature effectively reached by 
the subsystem, by fitting the Boltzmann distribution oc to the asymptotic 

{pniEn)}n=o,...io values. For this particular example, one finds that Tg^h = 0.99 ~ Tbath- 


White noise - Harmonic Fext 


Pn (t—>oo) Tsuh 



Fbath 


Figure 8: The asymptotic distribution of 
the eigenstate weights Pn=o,...io (red dots), 
obtained with A = 0.5 and Tbath = 1; 
function of the corresponding eigenenergies 
A„=o,...io- It fits the Boltzmann distribu¬ 
tion (oc 6“^/^="*=) with Tsnh = 0.99 (dashed 
line). 


Figure 9: Asymptotic subsystem tempera¬ 
ture Tsub as a function of the bath temper¬ 
ature Tbath for two different friction coeffi¬ 
cients: A = 0.1 (red circles) and A = 1.5 
(blue crosses) corresponding respectively to 
a weak and strong coupling. The dashed 
line corresponds to the ideal case Tgub = 
Tbath ■ 
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In Fig. we compare the temperature reached by our subsystem Tgub to the bath 
temperature Tbath (defined as the temperature entering the noise correlation). For a large 
range of temperatures and independently of the friction coefficient A and initial state, 
we observe that Tg^b — ^bath and that the asymptotic distributions of the weights are 
Boltzmannian. One can thus conclude that the subsystem correctly thermalizes when one 


uses the white noise autocorrelation (13). This is consistent with the assumption made 


in section 4.1 as regards the asymptotic behavior of any initial state and the stochastic 
evolution of its centroids. In this respect, Figs, [^and [^can also be understood as sanity 
tests of our numerical procedures. 

The total uncertainty on the asymptotic values, for a statistic of a few thousands 
of realizations, grows with the temperature from ~ 2% at Tbath = 0.1 to ~ 10% at 
Tbath = 5. An additional averaging over a time range Af once the equilibrium reached, 
leads to more reliable results whose accuracy then follows the common statistical law 
OC l/V^T-stat X At'. 


To conclude, we have generalized the analytic results obtained in Sec. 4.3.1 with an 
initial Gaussian wave packet to other initial states. We can thus conjecture that within 


the case of the harmonic potential and the white noise autocorrelation (13), the SLE 


universally leads to the thermal equilibrium of statistical mechanics. By universal, we 
mean that this result is independent of the bath temperature, friction coefficient and 
the initial state. Although this result was only expected at the weak coupling limit (as 
explained in the introduction), it is also reached in the intermediate and strong coupling 
regimes. 

4-4- Equilibration with a colored noise 


We now use the colored noise (15) to study the SLE ability to bring a subsystem 


to the thermal equilibrium of statistical mechanics. We concentrate our study on the 
asymptotic distributions of weights, as other properties are similar to the ones discussed 
in the previous section devoted to the white noise case. Eor initial Gaussian wave packets. 


the equations governing the evolution of a, Xc\ and pci established in section |4.1| hold 
as well, so that each realization of the evolution with a stochastic force leads to an 
asymptotic wave packet of width s/hJmujQ. It is thus possible to extend the method 


used in section 4.3.1 Eor this purpose, one needs to evaluate the asymptotic distribution 
W (xci,Pci) for the case of a classical damped harmonic oscillator driven by some colored 
noise. This can be done in several steps. One first decouples the equations of motion 


(18) by setting 


Xcl(t) 

Pc\{t) 


= c+(t)u+ + C-{t)v- 


where and V- are eigenvectors of the matrix 

1 

—mX± 


0 —1/m 


mUJn 


A 


(36) 


with respective eigen¬ 


values A± = -2 ± V - CUq. v± = 


The equations of motion for c± then 
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write 


c+ — -A+C+ - 


Fr 


and c- = —X-C- + 


Fr 


my/ A'^ — 4 (jJq my/ 


(37) 


Each of these equations describes the motion of a free particle undergoing ohmic friction 
and colored noise and admits a trivial solution, e.g. 


c+{t) = - 


1 


m^A? — 4 u!q Jo 




FR{t') 


(38) 


In the asymptotic limit, the autocorrelations between Xc\ and pd can thus be generated 
from those between c+ and c_. In particular, one obtains 

Sd) ) = (4)^+®^+ + (c+C-)(u^®u_ + u^®u+) + (c?.)u^®u_. (39) 


{XclPd) 

Using equation (38) as well as the corresponding one for c_ evolution, one obtains 

r*+oo /*H-oo 


(4) 


1 


asympt 


m?{A^ — 4 cc;q) Jq 


dti 


dt2 e"^±(‘i+*2)C(ti - t2 ), 


(40) 


10 


where C'(r) is the force autocorrelation {FR{t)FR{t + t)). Proceeding to the variable 
changes S = U +12 and t = |ti — t 2 | leads to 


(c^) asympt — 


71c(A±) 


m‘^{A‘^ — 4a;Q)A: 


di 


(41) 


where Cc designates the Laplace transform of C, while similar calculation for {xc\Pc\) 
gives 

'^c(A+) + Cc{X-) 


(XdPcl)asympt - ' ^2^(^2 _ 4 ^ 2 ) 

Gathering all results, the asymptotic correlations write 

1 


(42) 


asympt — 


(A+£c(A-)-A_£c(A+)) 


Am^ojQ y/ — 4 loq 
(Pcl)asympt = Ay/ A 4 ^ 4+'^c(A+) — X-Fc(X-)) 

(^dPd) asympt — 0. 


(43) 


For A < 2 ljo one obtains complex eigenvalues with however A_ = A+ and £c(A_) = 
jCc(X+), which guarantees that both (x^[)asympt and (p^i)asympt are genuine real quanti¬ 
ties: 


(3:cl)asympt — 


(Pcl) asympt — 


1 f 3f?£c(A+) aEc(A+) 


m^uj, 


2,.,2 




x/^^i - 41^ 


3f?£c(A+) , 9/:c(A+) 






(44) 
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In the white noise case limit a —?• 0, one has C'(r) = B5{t) with B = 2m^A;Tbath ac¬ 
cording to Einstein relatiorj^ As a consequence, Cc{^±) = w-^^Tbath and one recovers 
(a;ci)asympt = as Well as (p^i)asympt = m/cTbath, SO that each degree of freedom 

carries the expected average energy For the general case we express 


/t2 \ ^ _ ^^bath 

\‘^cl/3'Sympt — o • X 

and (Pcl)3'Sympt — '^^^^bath ^ •> 

mwg 

where 


A+£c(A-)-A_£c(A+) 

„ A+£c(A+) - A_£c(A-) 

mA/cTbathV^^ - 4a;^ 

" mAkTo.,^^A^ - Aul 


(45) 


(46) 


are reduction factors which encodes the deviations with respect to the white noise 
case. Finally, one easily demonstrates that only the second cumulants (x^j)asympt and 
(3;ci)asympt are non-zero for stochastic forces Fr of Gaussian nature. The ensuing asymp¬ 
totic distribution of the centroids (xci,Pci) thus admits the simple form 


2 2 2 

_ Pci 

W {Xc\iPc\) OC e 2mfcTi,ath’’P . 


(47) 


For the noise autocorrelation © , one has 


^c(A) 


mA f°° huo ( ^ 

vr Jq ex.-p{hw/kT^iMi) - 1 


oj — iX 


did 


(48) 


so that rx and only depend on the dimensionless quantities AiTbath/hwo and A/ujq. 
To our knowledge, the integral (48) has no simple expression. It however admits two 
interesting limiting cases: 


• High temperature: When feTbath 3> h|A±|, exp(^a;/A:Tbath) can be approximated 
by 1 -|- hw/ZcTbath in the integral, leading to Cc{X) k. mAfeTbath and ss 1: 

one reaches the classical limit. 


• Brownian hierarchy: When such hierarchy is satisfied A~^ should be larger than 
any time scale of the subsystem. Under such conditions, A <C wq, A+ ss iojq + ^ 
and 


^c(A+) 


mA huj / i 

vr Jo ex.p{hj/kToa.th) - 1 + iA/2 


i 

UJ +OJo- 



doj. 


(49) 

One has A <C AiTbath/h as well and the poles ±(ti;o — iAl2) can then be considered 
to be infinitely close to the real axis, with — *P-v. ^ 


One easily checks that this is true for the colored noise (151 when Tbath —> -too. 
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and 


— %A.j^ 


= resulting in 


^c{^+) 


mAhujo 


+ 


iraA 


exp(;ia;o/A;Tbath) - 1 ' tt Jq exp(;iw//crbath) - 1 - wg 


f 


hui 


2uJo 


duj. 
(50) 

According to equation ( |44[ ), the small A factor showing up in the real part of 
Cc{^+) is compensated by an equivalent factor in the denominator. Therefore 

{pV 


(2^cl)asympt — 


M. asympt 

2,,,2 


m^LJ, 


0 


and 


1 hujQ 

mujl exp(;ia;o/fcrbath) - 1 

^o/^i?bath 


exp(;ia;o/A:rbath) - 1 ’ 


(51) 

(52) 


On Fig. 10 we illustrate the quantities Vx and Vp as a function of feTbath/^o for sev¬ 
eral values of AjojQ. One notices that < 1 for finite A/ojq, which breaks the 

equipartition of energy. 

Tr & Tn 



'bath 


Figure 10: The and Vp reduction factors for several friction coefficients A (in units of wq), as a 
function of Tbath (in units of fujo/k). One has generically rx{A',T) > rx{A,T) and rp(A',T) < 
rp{A,T) for A' > A, as represented by the arrows. 


Following the method of section 4.3.1, the asymptotic weight of the state writes 

(_l)n Qn 


Pn = 


n\ drj 

where the generating function G is defined as 


G{u)q, a, Tbath) 


(53) 


r?=l 


^(wo, A,Tbath) = / lF(xci,Pci)e 




2 i j j _ (ii)i^:p (ii) 

dxcidpci — j j 


nj {r])U^ (t]) 


with 


^xid) =V + and Up{r]) =7]+ 


(54) 


(55) 
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Quite generally, and the pn generated from G are more involved than the 

simple power law oc found in equation (32) for the white noise case. As a 
consequence, deviations from usual Boltzmann distributions are expected for the pn- 
Noticeable exceptions are: 


• the case of a classical noise obtained for large T^athi for which = 1 and 

Pn oc exp (^-nln (^1 + ^1^)) ~ exp (-^In/A^Tbath) 


the weak coupling case obtained at small A for which, according to equation (52), 
IIj, = lip = r/ - 1 + exp(;i(uo/^7bath) implying that oc exp(-n/iwo/^7bath) oc 
exp(-£;„/A;rbath) as well. 


In these two cases, the temperature governing the distribution of weights for the subsys¬ 
tem (Tsab) is found to be exactly the bath temperature Tbath- 


Pnii^co) 7?„(t^Oo) 



Figure 11: The asymptotic distributions of the eigenstate weights Pn=o,...8 (joined by lines) 
function of the eigenenergies A„=o,...8) obtained with different friction coefficients, measured in 
units of ojo- A = 0.1 (solid lines), A = 0.5 (dashed lines) and A = 1.5 (dot-dashed lines) and 
temperatures (in units of huJo/k): Tbath = 0.2 (left), Tbath = 0.5 and 1 (right). They are compared 
to the corresponding “ideal” Boltzmann distributions oc (thin lines). 


Those various aspects are well illustrated on Fig. [m where deviations with respect 
to the Boltzmann distribution oc exp(—i?a/Tbath) are observed for ’’large” friction 
coefficients or ’’small” temperatures. For these cases, the distribution of weights appears 
to be better described by the law pn oc exp(—T^/Tsub), i-e. by introducing some effective 
temperature Tguh specific to the subsystem, as explained in the introduction. Yet, the 
definition of Tg^h is not unique as the {pn} shows genuine deviations from a power law. 
As one is often interested in the low lying eigenstates in phenomenology (the fundamental 
and few lower excited eigenstates), a bona fide choice for Tguh will be adopted here as 


Tsuh {{Pn}) = 


El — Eq 

ln(pi/po) ’ 


(56) 
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For the harmonic potential, po and pi are evaluated thanks to equation (53) and defini- 

1 


tion (56) leads to 


Tsuh = -- 


ln( 1 


1 


1 + 


+ 


1 


^bath’’a: 


1 + 


^bath’'P 


(57) 



Figure 12: Asymptotic subsystem temperature Tsuh as a function of the bath temperature Tbath 
- both in units of huio/k - for three different friction coefficients (in units of loq): A = 0.1 
(solid line), A = 0.5 (dashed line) and A = 1.5 (dot-dashed line) corresponding respectively to a 
weak, intermediate and strong coupling. The straight dashed line corresponds to the ideal case 
Tsub — Tbath ■ 


In Fig. 12, we compare this effective temperature actually reached by our subsystem 
Tsub to the bath temperature Tbath used as input of the noise. When Tbath ^0.5, the 
subsystem thermalizes with Tguh ~ Tbath in a good approximation (we note a slight 
dependence vs the friction coefficient: the larger A the smaller Tsuh)- At lower tempera¬ 
tures, some discrepancies appear in the form of a saturation of Tsuh- If A < wq the onset 
of those deviations is however delayed until Tbath ^ A where the Brownian hierarchy is 


then broken. 

Using our numerical tool, we have also studied the case where the initial state is 
not a Gaussian wave packet. Our conclusion is that the asymptotic states always turn 
out to be Gaussian wave packets of width y^h/mwo (which corresponds to the ground 
state width). We thus conjecture that the main results obtained in this section for the 
harmonic potential are independent of the initial state. 


5. Equilibration with a linear oscillator 

In this section, we study the thermal relaxation given by the SLE ([^ with the linear 
ID potential Uext = Ki \x\/‘^ with the white (13) or colored noise. The linear 
oscillator allows us to test the SLE for a non-harmonic situation. The analytic resolution 
being far from obvious in this case, we perform the analyses resorting to numerical 
simulations. For all initial states investigated, we have found similar asymptotic features 
and we will therefore not discuss the role of the initial state any further. 
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5.1. Equilibration with a white noise 


As explained in Sec. the white noise (13) was initially derived for a harmonic 
potential. In this section, we test its ability to be extended to other types of potentials 
through the example of the linear potential. In the white noise expression (13), we 
set Eq to the corresponding ground state energy, i.e. Eq 0.509 hwo, with huo = 


(Ho) 



Figure 13: Numerical average energy {Hq) evolutions for different friction coefficients A (solid 
curves) and the theoretical evolution given by (20) with Tbath = 1, {Ho){t —>■ oo) = 1.52 and 
A = 0.1 (dashed curve). 


We first consider the average energy observable (Hq). As shown in Fig. 13 
asymptotic value of (Hq) differs in general from the expected value 


{Ho){t oo) = 


E^E^e 


--Fi/rbath 




1.52, 


the 


(58) 


and exhibits a strong A-dependence, in contrast with the harmonic potential case. For 
small friction coefficients (A < 1), the average energy evolutions are nevertheless in good 
agreement with the exponential rate (20) when one takes the actual (i7o)(t —>• oo) and 
effective Agg ~ A/2 values. 


Independent of the initial state, the asymptotic distributions of the weights Pn=o,...io 
are close to the Boltzmann distributions oc only when 1 < Tbath ^ 2 at weak 

couplings (see Fig. 14). At low temperatures strong discrepancies are observed: the 
higher excited states exceed the Boltzmann law, exhibit an alternating pattern and 
saturate at low weights. A dependence on the friction coefficient value is observed from 
the 2^^*^ (4*^) excited state at low (medium) temperatures. The latter explains the {Hq) 
dependence on the friction coefficient observed in Fig. [T^ a smaller friction coefficient 
is observed to generate higher populations for the excited eigenstates and thus a higher 
average energy. At large temperatures Tbath ^ 5, the distribution of weights is difficult 
to evaluate because of statistical fluctuations and numerical scheme imperfections. 


22 

















Pn (t^oo) p„ (t^oo) 




Figure 14: The asymptotic distributions of the eigenstate weights p„=o,...io (joined by lines) 
function of the eigenenergies i5n=o,...io (vertical lines), obtained with different friction coefficients 
A = 0.1 (solid lines), A = 0.5 (dashed lines) and A = 1.5 (dot-dashed lines) and temperatures 
Tbath = 0.2 (left), Tbath = 0.5 and 1 (right). They are compared to the corresponding “ideal” 
Boltzmann distributions cx 6“^/^*=“*** (thin lines). 


Given that the SLE does not lead to genuine Boltzmann distributions in this case, 
we will resort our definition (56) of the subsystem effective temperature which is 
obtained by fitting the Boltzmann law to the first two weights (pq and pi) in lin-log 
space. The evaluation of such Tgub vs. Tbath is shown in Fig. [T^ For any value of the 
friction coefficient, one observes clear deviations from the “ideal” Tgub = Tbath line at 
low and high temperatures, which could be interpreted as some inefficiency of the SLE 
to heat up the subsystem. In particular, one does not recover this identity even in the 
small coupling limit. It appears in fact that the Tgub(Tbath) law is rather insensitive to 
A. 


White noise - Linear 14x1 

Tub 



Figure 15: Asymptotic subsystem temperature Tguh as a function of the bath temperature 
Tbath for two different friction coefficients A = 0.1 (red circles) and A = 1.5 (blue crosses) 
corresponding respectively to weak and strong couplings. The dashed line corresponds to the 
ideal case Tgub = Tbath- At high temperatures our accuracy on Tguh is low due to a very large 


time required to reach the asymptotic behavior and a large uncertainty as in Sec. 4.3.2 
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In Tab.[^ we assess the relevance of the effective Boltzmann laws obtained 

once Tsub has been defined, by counting the numbers of low lying eigenstates which are 
close enough to this law. At low temperature, only a limited set of weights are found to 
be encompassed by this law. 


— Number of weights close to oc e ^/Tub ? — 

Tbath \ A 

Small 

Intermediate 

Large 

Low (Tbath < 0.5) 

3 

2 

2 

Medium (0.5 < Tbath < 2) 

5 

5 

4 

High (Tbath > 2) 

10 

9 

8 


Table 1: Approximate number of weights close to the corresponding Boltzmannian oc 

One can consider the agreement to be poor from 2 to 4 weights, good from 5 to 7 and very good 

from 8 to 11. A better agreement is obtained toward the weak and/or high temperature regimes. 


In view of these elements, we conclude that the white noise (13) is not quite suitable 
to obtain an acceptable thermal equilibrium (in the sense of pn oc with other 

potentials than the harmonic one. Nevertheless, if one is interested in a limited number 
of low lying eigenstates (see Tab. [^, this formalism could be used for phenomenological 
purposes by performing a rescaling in the noise expression (13): either by changing the 
value of Eq (optimal value is found to be Eq = 0.33) or by choosing the input Tbath such 
as to obtain the desired Tsub = Tbath- Conversely, this study confirms the very specihc 
nature of the harmonic potential upon which general conclusions should not be drawn 
as regards the applicability of any scheme aiming at describing the thermalization of 
quantum subsystem. 


5.2. Equilibration with a colored noise 

Unlike the white noise, the colored noise (15) was derived without assumptions on the 
potential. In this section, we test its ability to be extended to other potentials through 
the example of the linear potential. 

As shown in Fig. 


16, the {Hq) average energies show similar features to the ones 


obtained with the white noise (Fig. 13), although one observes some overshooting of the 
expected value (1.52) at small coupling. The asymptotic distributions of the weights 
are found to be independent of the initial state and are displayed in Fig. [m They are 
given with the same values of {A, Tbath} for the harmonic oscillator case (Fig. 11) 
which constitutes the reference point for our analysis. In that case, one has found two 
overlapping regimes for which the Boltzmann distribution is recovered, namely the high 
temperature regime Tbath 1 and the weak coupling regime A <C {1, Tbath}- For the 
linear potential, one observes as well a good agreement with the Boltzmann distribution 
in the high temperature regime, as illustrated by the Tbath = 1 set of curves in the right 
panel of Fig. 17 with larger dispersion than for the harmonic potential, though. 

For a fixed Tbath ^ 1, large deviations are always found in the strong coupling regime 
A > 1 , especially when Tbath 1 (as illustrated by the dashed curve on the left panel). 
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Figure 16: Numerical average energy {Hq) evolutions for different friction coefficients A (solid 
curves) and the theoretical evolution given by (20) with Tbath = 1, {Ho){t —> oo) = 1.52, and 
A = Q.l (dashed curve). 


Pn (t^oo) p„ (t^oo) 




Figure 17: The asymptotic distributions of the eigenstate weights Pn=o,...io (joined by lines) 
function of the eigenenergies i?n=o,...io (vertical lines), obtained with different friction coefficients 
A = 0.1 (solid lines), A = 0.5 (dashed lines) and A = 1.5 (dot-dashed lines) and temperatures 
Tbath = 0.2 (left), Tbath = 0.5 and 1 {right). They are compared to the corresponding “ideal” 
Boltzmann distributions (x 6“^/^*=“*** (thin lines). 


In the weak coupling regime, however, the Boltzmann distributions are only matched for 
a limited number of low lying eigenstates which appears to be reduced as compared to 
the harmonic case. It is still an open question whether one has lim^_^oTn oc 
in the weak coupling limit (a convergence that could be at best non uniform) or if 
finite deviations from the Boltzmann distribution survive. To conclude our comparison, 
one should note that the distributions of weights for the linear potential show overall 
richer pattern at small Tpath, as for instance non monotonous dependences vs the friction 
coefficient or similar alternating patterns as in the white noise case (see Fig. 14) with 
however lighter oscillations. 

Despite these discrepancies, the relation Ts^h vs. Tbath (Fig. 18), obtained by focusing 
on the lowest excited states, is interestingly close to the one obtained with the harmonic 
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potential. In the weak coupling regime, one recovers in particular ~ ^bath) except 
for Tbath ^ 2 where one naturally recovers the white noise results, with Tguh > Tbathj as 


already observed on Fig. 15) . These observations confirm the rather general nature of 


the colored noise (15), which might thus be combined with a wider class of potentials and 


used in a good approximation for thermalization studies especially in the weak coupling 


case. 


Colored noise - Linear Fext 

Tub 



Figure 18: Asymptotic subsystem temperature Tguh as a function of the bath temperature Tbath 
for three different friction coefficients: A = 0.05 (red circles), A — 0.5 (green diamonds) and 
A = 1.5 (blue crosses) corresponding respectively to a weak, intermediate and strong coupling. 
The dashed line corresponds to the ideal case Tgub = Tbath- 


5.3. Application to an interesting subatomic system 

Although one has strong evidence that protons, neutrons and other hadrons are made 
of more fundamental objects - the quarks and the gluons - these constituents are usually 
confined inside those hadrons. At high temperature or density, the coupling constant 
is reduced and the quarks and gluons are expected to be found in a deconfined phase 
(the so called quark-gluon plasma (QGP)), probably achieved for a short lapse of time 
in the early universe. Nowadays, it is thought that such a state could be formed and 
investigated at best in ultra relativistic heavy ion collisions (URHIC) taking place in 
large ion colliders, as for instance RHIC and LHC located respectively at Brookhaven 
National Laboratory and at CERN. One of the possible signatures of the QGP is the 
melting of all hadrons made of cc and bb quark-antiquark pairs - the so-called quarkonia 
otherwise quite robust in less extreme conditions - resulting in a suppression of their 
production as compared to a situation were no QGP would have been formed. Such a 
mechanism was postulated in 1986 by Matsui and Satz m and was indeed observed 
later on, but the analysis is however not quite conclusive (see |73! for a recent review), 
one of the main reason, on the theory side, being that the dynamics of those quantum 
object is not properly taken into account in most models. It is the purpose of this section 
to investigate whether the SLE could be appropriate to deal with such a system. 

Quarkonia can be described in non-relativistic quantum mechanics through phe¬ 
nomenological potentials, the most simple of them being the Cornell potential (see m 
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and references therein): 


Vcomell(?’) :=-h Sr , 


(59) 


where r is the distance between the quark and the antiquark, k = 0.52hc and s ~ 1 
GeV/fm is the string tension. It is thus the relative coordinates which should be taken 
as the fundamental degrees of freedom of the SLE. To make the link with the simple 
ID linear potential studied in sections 5.1 and 5.2, one can note that the Coulomb part 
oc in (59) mainly acts on the ground state but has much smaller influence on excited 
states, especially for cc bound states (charmonia). For the sake of the analysis, we will 
thus associate Vcorneii(’’) to a linear ID potential with Ki = 2s ^ 2 GeV/fm. As for 
the quark mass, the optimal parameters extracted from spectroscopy analysis m are 
nic ~ 1.84 GeV/c^ and rrih ~ 5.2 GeV/c^. Concentrating on the charmonia case, one 
gets a reduced mass /r Ri 0.9 GeV/c^ leading to 


hwQ = 


{hcKiY 

fJ.C^ 


(0.2 X 2)2 
041 


0.55 GeV, 


(60) 


as he ~ 0.2 GeV • fm. Typical temperatures reached in most energetic URHIC are of 
the order of /cTbath £ [0-15 GeV, 0.6 GeV], resulting in G [0.3,1.1]. 

As for the friction coefficient of heavy quarks interacting with a quark gluon plasma, 
one can find several estimations in the literature, and we will here refer to the calculation 
of one of the authors m- In this work, it was found that A ~ 1.5A:rbath(c/fm), with 
^^bath is expressed in GeV and where the increase with the temperature results from the 
increasing density of interacting particles (quarks and gluons) which are the ground for 
friction. Consequently, one finds HA ps 1.5A:rbath|^ ~ O.SfcTbathj implying that one never 
encounters a breakdown of the Brownian hierarchy, even at the smallest temperatures 
achieved in URHIC. As HA < /cTbath ^ lAhcoo, we infer that this system lies in the 


intermediate coupling regime. Examination of Fig. 17 reveals that states up to the 3rd 
excited one {E 3 2huJo) should be correctly described by the SLE, which represents most 
of the charmonia spectroscopy. All together, we conclude that the SLE is a relevant 
approach to study the dynamics of quarkonia formation in URHIC and probe more 
accurately the QGP formed in these collisions. 


6. Discussion and conclusion 

For the purpose of finding an effective formalism suitable to phenomenological appli¬ 
cations of open quantum systems, we have focused on the Schrodinger-Langevin equation 
Q. Its nonlinear friction term is commonly believed to maintain the stationarity of the 
excited states of the uncoupled Hamiltonian Hq. We have shown in Sec. that the 
Madelung/polar transformation of the wave function leads to a nonzero damping for 
these states. In this way, we have reconciled the SLE with the intuitive expectation that 
the dissipation process should act on any state in order to bring the subsystem to its 
ground state. 
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We have then focused on the question of the thermal relaxation dynamics given by 
the SLE for two different potentials and with two different noise operators, taken as c- 
numbers: the white noise ( 13 ) - which has been derived by Senitzky m and subtracted 
by its term of ground state fluctuations - and the colored noise (15) derived by Ford, 
Kac and Mazur m- 

We first considered the case of a harmonic potential for which most of the results 
regarding asymptotic states and distributions could be established analytically for both 
kind of noises. When the subsystem undergoes a white noise, the SLE has demonstrated 
its ability to bring any initial state to the thermal equilibrium of statistical mechanics 
(i.e. Boltzmann distributions of the uncoupled subsystem energy states), irrespective of 
the coupling strength and of the bath temperature - confirming the assumption made 
by Messer [Sjj - although such equilibration is generally expected only in the weak cou¬ 
pling limit (as explained in the introduction). For a colored noise, exact thermalization 
toward the Boltzmann distribution was established only in the large temperature regime 
(fc^bath ^wo) or in the so-called weak coupling regime {A <C {woi ^^bath/^}) where 
the Brownian hierarchy is satisfied. Outside of these regimes, it has been shown that 
most of the deviations with respect to the Boltzmann distribution exp(—E„//cTbath) 
could be accommodated by introducing some effective temperature T^nh corresponding 
to the subsystem internal equilibration. The disagreements between Tgub and the bath 
temperature Tbath (input of the noise) observed at low temperatures or for large friction 
coefficient A can be attributed to the breaking of the Brownian hierarchy. 

The study for the case of a subsystem submitted to a linear potential was performed 
thanks to numerical investigations. While a rather clean equilibration was found in the 
high temperature regime, genuine non-Boltzmannian behaviors for higher excited states 
(which cannot be accommodated by a simple change of the temperature Tbath —>■ T^nh) 
and stronger friction coefficient dependences have been observed at low and medium 
temperatures for both kind of noises. Concentrating on a smaller subset of low ly¬ 
ing eigenstates, the colored noise has nevertheless led to better results in the sense of 
statistical mechanics, confirming its more universal nature. In particular, the effective 
temperature Tgub was found to follow similar behavior as the one introduced in the har¬ 
monic potential case and to converge toward Tbath in tke weak coupling regime. We 
thus conclude that the colored noise should be used preferentially, especially at low 
temperatures (under the condition that the Brownian hierarchy is preserved). 

The SLE and the quasiclassical Langevin equation seem therefore to share a common 
difficulty in the description of dissipative evolutions outside the nearly harmonic and 
free potential cases [TTj and especially outside the classical high temperature limit [63j . 
Nevertheless, further analyses - such as in quantum tunneling [69| - would be required 
to establish a common behavior. 

If one focuses on applications where only the lower states are considered, a phe¬ 
nomenological model could be to rectify the observed differences between Tgub and Tbath 
by choosing an effective heat-bath temperature Tbath such as to reach the desired subsys¬ 
tem temperature Tgub = Tbath- It just requires the proper knowledge of the Tgub (Tbath) 
function as displayed in Figs. and Though one has to adapt this rescaling 
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to each situation, the SLE can thus be considered, in a good approximation, as a pos¬ 
sible effective alternative to quantum master equations and SSE for phenomenology of 
complex quantum systems (such as ions quantum transport, thermalization of quarkonia 
in a quark-gluon plasma, of nucleons in a nucleus,...Dealing with the full hierarchy 
of states in the general case would possibly require either the use of the colored noise 


correlation (14), of a q-number noise operator or of a more refined quantum treatment of 


the subsystem interactions with the heat bath. Establishing an Einstein relation specific 
to the SLE would also be an interesting perspective for this field. 

It should be noted that our analysis relies on the hypothesis that the asymptotic 
distribution of subsystem-eigenstates weights pn must be Boltzmannian whatever the 
potential and the coupling strength to the rest of the system (the heat bath). To our 
knowledge, such an assumption has not been universally established from fundamental 
principles (i.e. starting from the distribution of the full-system eigenstates and tracing 
out the heat-bath degrees of freedom) and should be considered more thoroughly in 
future studies as it could partially alter our conclusions. 


Appendix A. Numerical generation of noises 

We describe here a numerical method used to sample numerically colored noises 


characterized by finite autocorrelation times such as (15). In a first step, one generates 
stochastic stationary Gaussian random variables fj with zero average and correlation 
{rjfj') = At5jj', where At is the time step of the numerical scheme. Then, one builds 
the Gaussian random force E at a time ti - and assumed to be constant over the time 
step [U, ti + At] - from the weighted sum 


+ 00 


Fi= Y, 


(A.l) 


]=-oo 


where the weights Wi-j depend only on the difference i — j to guarantee the stationarity 
of the process. Then, the average of Fi is null and its covariance is given by 


Fi Fi,) = 


+ 00 


H-oo 


'Y Wi-jWi,-j,{fjrj,)= Yj Wi-jWi,-jAt, 


(A.2) 


] = -0O 


which, in the continuous limit At - 
{FR{t)FR{t')) = 


0, becomes 

r+oo 

W{t - t") W{t' - t") dt” , 


(A.3) 


with Wi = W(ti) for a given time step At. Then, one easily shows that the Eourier 
transform of W is just the square root of the power spectrum P{uj) of the retained 
noises, i.e. 


Piuj) = |W(a;)|2 


(A.4) 
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with 


P{ijj) = 2 mA 


hoj 


for the colored quantum noise (15) and 


exp(;ia;/A:rbath) - 1 ’ 


P{ijj) = lim B exp ( — 


O'—^-O 


(A.5) 


(A.6) 


for a Gaussian C with autocorrelation time a. the white quantum noise (13) is obtained 
from the latter when n —0, but in practice it is sufficient to take a much smaller then 
the typical times governing the subsystem evolution. Then, one gets explicitly 


1 

W(r) = - / ) cos (cur) dw 

^ JQ 


(A.7) 


and then the stochastic variables {T)} through equation (A.l). 


pR 



{Fr(\)FAX+t)) 



Figure A. 19: Left: Example of one colored noise (15) realization obtained with the described 
numerical method. Right: Corresponding analytical (dashed black curve) vs. numerical (orange 
dots) covariances over time. 


In Fig. A. 19 (left), an example of a colored noise realization obtained with the 
described numerical method is shown. In Fig. A.19| (right), the corresponding numerical 
correlation over time is successfully compared to the analytical expectation. 

Besides, one can easily show that the variables defined in this way are Gaussian. 
Similar algorithms can be found in the literature and have been successfully used in SSE 
and other formalisms (see I57j and references therein). 


Appendix B. Evolution of a Gaussian wave packet under the SLE with a 
harmonic potential 

The evolution of a Gaussian wave packet under the SLE has been extensively dis¬ 
cussed by several authors (see [32] and references therein). Here we derive anew the 
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essential results for our study, concentrating on the asymptotic shape of the solution. 
We consider the ID SLE with a harmonic potential: 

^ h^(5'(a:) - (5)) V' - FR{t)x'il) + (B.l) 

where K = mwg, A is expressed in s“^ and plays the role of a relaxation rate, while S 
is the phase of V' and FR{t) is a stochastic force. We consider the following Ansatz as a 
specific class of solutions 


oc exp 


i 

h 


{a{t){x - Xci{t)f +pciit){x 


Xc\{t)) + lit)) , 


(B.2) 


where a governs the wave packet width (with Q(a) > 0), Xd is the centroid in usual 
space and pd is the centroid in momentum space, both assumed to be real. 7 combines 
a phase factor (real part) and an normalizing factor (imaginary part). The l.h.s. of the 
SLE gives 

= - [a{t){x - Xci{t) f+ 2a{t)xci{t){x - Xci{t))+pci{t){x - Xci{t))- 

Pci{t)xci{t) + 7(t)] ipA ■ (B.3) 

As for the r.h.s, one has 

” ^ ^ {-2iha(t) + [2a{t){x - Xd{t)) +Pci{t)f} (B.4) 

while the friction term leads to 

hA{S{x) - (S)) ipA = A [3f?(a(t))(x - Xci(t))^ +Pci(i)(ic - Xci{t))] tpA ■ (B.5) 


Both l.h.s. and r.h.s. correspond to pJA multiplied by second degree polynomials. Equat¬ 
ing terms oc leads to a first equation: 



A3?(a) 


mcjQ 

2 


(B.6) 


Writing 

"^3:2 = [(x - Xci(t))2 -h 2xci(t)(x - Xcl(t)) -h Xd(t)] , (B.7) 


all terms oc 
terms oc (x 


(x — Xci(t))2 cancel provided equation (B.6) is satisfied, 

- Xci(t)): 






- [2a{t)xd{t) -bPci(0] (x 

/ 2Q;(t) 


-Xc\{t)) = 

+ A] pci{t) - Fnit) + mujQXci{t) 


(x 


Xcl{t)) , 


(B.8) 
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where we have used Fji{t)x = Ffi{t){x — Xci{t)) + FR{t)xci{t). As a{t) G C while other 
quantities are real, real solutions for Xd and pd are obtained by identifying the powers 
of a in both members of the equation, leading to 


Xci = — and Pd = -Apd + FR{t) - mulxd 
m 


(B.9) 


which are just the equations of motion for a classical particle evolving in some harmonic 
potential under the action of a friction term and a fluctuating force. Finally, we are left 
with the constant terms in the polynomial: 


Pc\{t)xd{t) - 7(t) = - 


iha{t) Pd{t) 


+ 


muj. 


m 


2 m 


- Fii{t)xd{t) + -rrXdit ), 


that is, using equation (B.9) 


i{t) = 


iha{t) pliit) 


+ 


muj, 


m 


2 m 


+ FR{t)xd{t) - —Xd{t). 


(B.IO) 


(B.ll) 


Expressing a = 51?(a)+i9‘(a) as well as 7 = 51?(7)+z9'(7) and equating real and imaginary 
quantities leads to 


d3f?(7) _ h9(a) Pd{t) 


dt 


m 


+ 


muj, 


and 


2 m 


+ FR{t)Xd{t) - T^x^iit) 


(B.12) 


dt m 

The latest equation represents nothing but the norm conservation. Indeed, one has 

29(7) 

[ dxiV'Ap = [ oc 


(B.13) 




and norm conservation reads 


2 d^{'y) 1 d^{a) 

h dt 23 ( 0 ;) dt 


= 0 , 


which can be recovered by combining the imaginary part of eq. (B. 6 ), i.e. 

43'(q;) 3f?(Q;) 


^(q;) = —- 


m 


(B.14) 


(B.15) 


(B.16) 


and equation (B.13). We now turn to the study of equation (B. 6 ), in particular its 
asymptotic behavior. It is obvious that this equation has a single stationary point, 
namely Ost = corresponding to the ground state width 7 ^- We therefore 

express a in units of {a = and t in units of {t = t/ujq) and obtain the 

reduced equation 

^ = —z^ — 1 — A^{z ), 

dr 
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where A = —. We first discuss the case ^ = 0 (no friction), for which the general 


solution of equation (B.13) writes. 


z(t) = 


z(0) COST — sinr 


(B.18) 


z{0) sinr — cost 

Analysis of this solution shows that z{t) revolves, in the complex plane, on a circle of 
center C = ^0, ^ and of radius 

V(1 + ^2(0))(i + ^2(0)) 


r(x,y) = 


29(z(0)) 


(B.19) 


which is an invariant of the motion. On any of these circles, one has max(S5(z)) x 
min(S5(2;)) = 1, which implies that the stationary solution Zgt = (0,i) is located inside 
each of them, but does not correspond to the asymptotic behavior of any non-trivial 
solution. We now turn to the effect of finite A on the evolution. Using the decomposition 


z{t) = x{t) + iy{t), the real and imaginary parts of the reduced equation (B.17) write 
= fix, y) with fix, y) = ( ^ ^ ) • (^-20) 


dr 


-2xy 


Even if this system admits no analytic solution to our knowledge, an illustration of / in 


the form of a vector field plot (see Fig. B.20) reveals its global rotating naturq^ around 
the center Zgt- 


y 



Figure B.20: Illustration of the vector field / defined in equation (B.201 for A = 1/4; the black disk 
corresponds to the stationary point Tst = (0,1). The spiral curve represents a solution of the equations 
of motion (|B.20[) for initial conditions (a:(0) = 0,j/(0) = 2). 


^^Conhrmed by rot/ = —Ay < 0 in the upper complex plane. 
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As compared to the A = 0 case for which the evolution driven by / takes place on a 
single circle (preserving the invariant r{x,y)), one can show that the —Ax terms causes 
some ’’damping” of the motion associated with a continuous decrease of 


dr"^ _ Ax^{1 + x"^ + y^) 
dr y‘^ 


This prevents z{t) to stay on any given circle. Instead, solutions for finite A are inward 


spirals which all ends at Zgt, as illustrated by the solid line in Fig. B.20 


lim z{t) = Zst ■ (B.22) 

r^+cxD 


This proves that all Gaussian wave packets evolving in some harmonic potential under 
the SLE ultimately acquire the ground-state width. 

To conclude, we wish to make the connection with the results established using the 
hydrodynamic form of the SLE (see for instance |32] for a recent discussion). In this 
approach a Gaussian Ansatz is made for the probability density 


1 


P = 


where a satisfies the Pinney equation 


exp 


{X - Xcl) 

2a{tY 


'd T Ad T ujqQj — 




(B.23) 


(B.24) 


for 


Comparing the Gaussian Ansatz (B.23) with our Ansatz for ijj, equation (B.2 ) lead s to 
^ ^ ~ d = Using, the equation of motion (B.16 

the imaginary part of a, one obtains 

8 A(a) 5?(Q;)a^ 2 3?(a)a 


a = 


hm 


m 


(B.25) 


In order to obtain a closed differential equation in a, we differentiate equation (B.25) 
and use the equation of motion (B.16) for the real part of a: 


d = — (3?(d)a + I?(a)d) = — 
m m 


— (A(a)^ — 5?(a)^) — A3?(a)-) a + I?(a)d 


m 


(B.26) 

We then perform the final substitutions I?(a) —)• ^ and 'A{a) —)■ ^ and indeed recover 
equation (B.24). 
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